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INTERACTION  BETWEEN  STEADY  NON-UNIFORM  TWO-DIMENSIONAL 
CURRENTS  AND  DIRECTIONAL  WIND-GENERATED  GRAVITY  WAVES 
WITH  APPLICATIONS  FOR  CURRENT  MEASUREMENTS 

I.  INTRODUCTION 

In  recent  years,  the  interactions  between  surface  waves 
and  variable  currents  have  been  studied  by  many  researchers. 
Unna  (1942)  and  Evans  (1955)  studied  this  phenomenon  without 
considering  the  effects  of  nonlinear  energy  transfer  between 
surface  waves  and  currents.  The  analysis  given  by  Drent 
(1959)  was  for  plane,  long-crested  waves  and  rectilinear 
shear  flow  in  the  absence  of  surface  tension  and  viscosity. 
Hughes  and  Stewart  (1961)  studied  the  same  case  by  experiments 
except  the  waves  are  in  the  gravity-capillary  range.  Longuet- 
Higgins  and  Stewart  (1961,  1964)  advanced  the  concept  of 
wave-current  interactions  by  introducing  nonlinear  energy 
transfer,  between  waves  and  current,  in  the  form  of  the 
inner  product  of  the  wave-induced  radiation  stresses  and 
gradients  of  the  current.  In  their  analysis  the  net  total 
energy  change  in  the  dynamic  system  was  assumed  to  be 
negligible. 

Waves  do  change  their  characteristics  both  kinematically 
and  dynamically  whenever  they  encounter  currents.  This 
phenomenon  of  changes  has  not  only  been  observed  by  mariners 
for  centuries  but  has  also  been  systematically  detected 
from  satellite  with  scanning  radiometer  as  reported  by 
Strong  and  DeRycke  (1973).  These  observations,  combined 
with  the  theoretical  treatise  on  wave-current  interactions, 
greatly  increase  the  possibility  of  inferring  currents  by 
using  remote  sensing  as  a means  to  observe  the  changes  in 
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ocean  wave  characteristics.  In  order  to  explore  this 
possibility  further,  Huang,  et  al.  (1972)  and  Tung  and  Huang 
(1973)  derived  quantitatively  the  measurable  physical 
parameters  such  as  dispersion  relationships,  wave  spectra, 
slope  spectra,  surface  roughness,  and  wave  height  distribu- 
tion for  wind-generated  gravity  waves  in  various  steady 
nonuniform  currents.  With  the  progress  made  recently  in  the 
field  of  remote  sensing,  these  physical  parameters  should  be 
readily  observable  from  moving  platforms  such  as  aircrafts 
and  satellites.  Furthermore,  based  upon  the  finding  of  Long 
and  Huang  (1976a,  1976b) , the  same  measurement  concept  and 
technique  can  be  extended  to  the  range  of  capillary-gravity 
waves.  Their  theoretical  predictions  were  confirmed  by  the 
laboratory  observations  on  the  changes  in  kinematic  and 
dynamic  wave  characteristics  for  this  spectral  range  of 
waves  with  variable  currents.  Their  studies  strengthen  the 
concept  of  determining  currents  by  observing  waves  remotely. 

However,  all  the  studies  mentioned  above  are  for  the 
cases  of  one-dimensional  waves.  Chen,  et  al.  (1973)  inves- 
tigated the  case  of  directional  fetch-limited  wind-generated 
gravity  waves  on  variable  shear  currents.  The  continuity 
equation  is  identically  satisfied  in  that  particular 
situation  where  the  current,  strictly  speaking,  is  still 
one-dimensional.  This  paper  will  discuss  the  case  where  the 
continuity  equation  is  not  necessarily  identically  satisfied 
and  will  provide  a complete  numerical  solution  of  the 
directional  spectrum  describing  the  interactions  of  wind- 
generated gravity  waves  with  two-dimensional  variable 
currents  in  deep  water.  The  method  may  be  extended  to 
include  more  general  cases.  Solutions  are  obtained  by  the 
method  of  characteristics,  reducing  the  partial  differential 
equations  to  ordinary  differential  equations,  the  solutions 
of  which  are  obtained  by  numerical  integration.  The 
directional  wave  spectra  generated  by  the  wind  and  measured 
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by  Longuet-Higgins , et  al.  (1963)  is  used  as  a guideline  for 
the  case  of  waves  with  no  current. 


II.  THEORY 

Let  a component  of  a general  wave  field  be  specified  by 
!(S,t)  = B(5,t)  (1) 


where  B(x,t)  is  the  amplitude  and  x(x,t)  the  phase  function; 
both  are  functions,  in  general,  of  position  vector,  x,  and 
time,  t-  The  wave-number  vector,  k,  and  wave  frequency,  n, 
can  be  defined  as 


k = 


n 


(2) 


From  the  two  expressions  in  Equation  (2)  it  implies  conse- 
quently that  the  wave  number  vector  is  irrotational  in 
space,  i.e. , 


v x jc  = o 


(3) 


5 0 (4) 

kinematic  conservation  law  of  waves.  The 
amplitude,  |k|,  is  defined  as  2tt/\  where  X is  the  wavelength. 
Wave  frequency,  n,  however,  is  not  the  same  as  in  an  intrinsic 
oscillatory  wave;  rather,  it  represents  the  total  wave 
frequency  which  is  actually  observed.  For  waves  fluctuating 
on  a current,  U(x,t),  this  total  wave  frequency  is  composed 
of  two  terms,  namely,  the  intrinsic  wave  frequency,  a. 


and 


3k  „ 

3t  + Vn 


which  is  the 
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i 
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and  the  convective  wave  frequency,  k • u(x,t),  as 

n = a + k • U (x,  t)  ( 

where  a and  jkj  satisfy  the  dispersion  relationship,  which 
is  defined  as 

a = g|k|  (( 

with  g the  gravitational  acceleration.  The  weak  nonlinear 
dynamic  interaction  is  neglected  in  Equation  (6) . 

From  Equation  (4)  it  follows  immediately  that,  for 
steady  state, 


which  implies  that  the  total  wave  frequency  is  an  invariant 
in  steady  state  condition.  With  Equations  (5)  and  (7),  we 
will  have 


I 


! 

I 


n - a + k • U(x,t)  = aQ  (8) 

where  the  subscript,  o,  indicates  the  case  when  there  is  no 
current.  From  the  kinematic  influences  shown  by  Equation  (8) 
the  "wave  frequency  dispersion"  will  produce  changes  in  wave 
characteristics . 

For  the  current,  which  is  averaged  over  the  depth,  the 
kinematic  conservation  equation, 


guarantees  the  condition  of  continuity  for  steady  state. 
Thus,  Equations  (3),  (8),  and  (9)  constitute  the  kinematic 
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conservation  conditions  for  the  combined  field  of  waves  with 
current. 

The  dynamic  influences  can  be  derived  and  evaluated  by 
the  energy  balance  equation  as 

^ ^ 3U  . 

V • E (U  + cg)  + ^ = £ i,j  = 1,2  (10) 

where  E is  the  wave  energy,  c is  the  group  velocity,  U.  is 
the  component  of  U in  x^-direction,  is  the  radiation 

stress  and  the  term,  e,  represents  the  "net"  energy  gained 
by  the  dynamic  system.  In  the  case  of  the  assumed  steady 
active  wind  field  the  energy  gained  by  waves  is  more  or  less 
equal  to  the  energy  lost  by  waves.  Thus  it  is  further 
assumed  here  that  £ is  negligibly  small  and  is  set  to  zero 
in  Equation  (10),  i.e. 

9U . 

V • E (U  + Sg)  ♦ S±j  =0  (11) 

This  equation  governs  the  energy  transfer  betweeen  waves  and 
steady  non-uniform  current. 

III.  NUMERICAL  SOLUTION 

The  coordinate  system  (x^,x2)  used  in  the  previous 
section  can  be  arbitrarily  chosen  so  that  the  coordinates 
can  better  describe  the  physical  problem.  In  other  words, 
the  terms  are  either  geometric  components  or  quantities.  In 
forming  the  numerical  solution,  the  Rectangular  Cartesian 
Coordinates  x^  and  x2  will  be  assumed  for  the  present,  as 
shown  in  Figure  1,  without  loss  of  generality.  Let  the 
x^-axis  be  the  direction  of  the  wind  velocity,  W,  and  let 
the  wave,  whose  wave  number  vector  is  k,  propagate  at  an 
angle,  9,  with  respect  to  the  x^-axis.  The  local  current 
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vector,  U,  is  making  an  angle,  <J) , with  wave  number  vector, 
k,  and  U is  also  moving  at  an  angle,  C,  with  respect  to  the 
x^-axis.  All  the  angles  are  positive  in  a counterclockwise 
direction.  With  the  configuration  shown  in  Figure  1,  the 
formulation  of  a numerical  solution  can  be  derived. 


A.  Kinematics 

Based  upon  Figure  1,  Equations  (3),  (8),  and  (9)  can  be 
rewritten  as 


sin9)  = ( | it  | cos9)  , 


(12) 


n = a + | k | | U | cosif  = cQ 


(13) 


and 


( 1 U | cos£) 


dq  <|S|  sin5> 


0 


(14) 


respectively.  A solution  of  Equations  (12)  and  (13)  deter 
mines  the  local  direction  of  propagation  with  angle,  8, 
for  wave  number  vector,  k.  Equation  (14)  describes  the 
characteristics  of  the  current,  U. 

Differentiating  Equation  (13)  with  respect  to  x^  and 
*2  t solving  for  99/9x^  and  96/3x2  in  terms  of  9]kj/3x^ 
and  3 |/3x2  and  combining  with  Equation  (12),  it  is  easy  to 

see  that  the  partial  differential  equation  in  the  form  of, 


F 


+ 


H 


(15) 


for  determining  k can  be  readily  obtained  where 
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F = cos9  - 


U 1 sin9  sintj) 


+ i U cosd) 


G = sine  + - ! °1  cos9  sin* 

I5!  cos< 


cos<(>  I sin0 


- COS0 


3 | U] 


+ | U | sin<£  ( cos0  I-2- 


35  ■ „ 35 

v-2—  - sin9  tt-2— 
3x2  ^X1 


(0.5|/g/Jk)  + JuJ  cos<p)  . (18) 

The  wave  number  ]k|  in  the  expressions  for  F,  G and  H is  a 
function  of  the  angles  9 and  $ and  is,  in  general,  a function 
of  the  space  coordinates  x^  and  X2*  Equation  (15)  is  solved 
by  the  method  of  characteristics  (see  Courant  and  Hilbert 
(1966)).  By  integrating  the  ordinary  differential  equations, 


= F, 


a family  of  one-parametered  characteristic  curves,  s,  can  be 
obtained.  With  subscript,  o,  denoting  the  initial  condition 
there  is  a unique  curve,  in  this  family  of  curves,  which 
passes  through  each  point,  (xq,  0Q) / on  the  initial  curve, 

C.  This  family  of  curves,  s,  generates  the  surface, 

9 = 9 (x,  xQ,  9 ) , which  is  the  solution  for  Equation  (15) 
with  initial  conditions. 

The  model  of  the  current  is  chosen  such  that  U is 
definable  in  the  positive  x^  half  space  which  includes  xQ 
or  (o, n) • To  obtain  the  local  direction  of  wave  number 
vector,  k,  at  a given  x.  Equation  (15)  is  integrated  by  the 
Runge-Kutta  numerical  method,  varying  the  quantity,  n,  in 
xQ  in  such  a way  so  that  the  characteristic  curve  passes 
through  the  given  x.  From  the  solution,  9 = 0 (x,  xq,  9 ) , 
the  magnitude  of  the  local  wave  number  vector,  |k|,  is 
determined  from  Equation  (13).  In  other  words,  in  order  to 
obtain  the  local  wave  number  vector,  k,  at  a particular 
space  point,  x,  an  initial  value,  n , for  xq  must  be  determined 
for  each  0Q  and  xQ  = (o,n) • Thus,  a set  of  n values  is 
obtained  for  a particular  space  point,  x.  This  set  of  n 
values  is  used  for  the  solution  of  the  energy  balance  equation. 

The  average  current  field  is  assumed  to  be  of  the  form 
| U 1 cos?  = a | X£ | , (22) 


and 


|U|  sin£  = bx^,  (23) 

where  a and  b are  constants  and  Equations  (22)  and  (23) 
satisfy  the  required  equation  of  continuity,  Equation  (14). 

A verification  of  the  numerical  formulation  is  made  by 
assuming  only  a transverse  gadient  for  the  current  velocity. 
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The  results  show  that  |k  sin9  = Ik  I sin9  for  each  x which 

11  'o'  o 

agrees  with  the  analytical  results  by  Longuet-Higgins  and 
Stewart  (1964). 


B.  Energy  Balance  Equation 


The  radiation  stress  tensor,  S^,  has  the  form  of 


S.  . = 
19 


( 


\ 


E 2a 
2 cos  9 


^ sin0  cos9 


E E 2 

sin9  cos0  sin  9 


(24) 


J 


for  waves  in  deep  water  by  referring  to  Figure  1.  Further- 
more, Equation  (11)  can  be  reduced  to 


3E_ 

3x. 


+ M 


3E_ 

3x, 


= N 


(25) 


where 


L = 4-jl  COS0  + I U I cos?  , 


(26; 


M 


->■ 

c 


j1-  sinQ  + I U | sin£  , 


(27; 


and 


N = - E 


1 / 3 c n i^i  . A 9 0 

2(3^  cose  - l°l  Sln0  sq 


+ I mi  cos?  - | U | sin?  35 


3x, 


3x, 


1 + 


2a 
cos  9 
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\ 


i 


*( 

( 


3 I c I . » i -►  t „ 36 

3^-1  sine  + |c|  cose  3J- 


il£J 

3x„ 


sin5  + I U I cos? 

y X 


) 


+ ^ sin9  cos0  ^ 1 cosC  - |u|  sin5 


+ fir-  sin?  + iffi cos?  M 


r)l 


(28) 


and  | c j is  the  phase  speed  of  the  wave  whose  wave  number  is 

Ikl. 


The  solution  for  Equation  (25)  is  obtained,  once  again, 
as  for  Equation  (15),  by  the  method  of  characteristics.  By 
integrating  the  ordinary  differential  equations. 


(29) 


dx2 

dp 


M, 


(30) 


and 


(31) 


the  family  of  one-parametered  characteristic  curves,  p,  can 
be  obtained.  There  is  one  and  only  one  of  this  family  of 
curves,  p,  which  will  pass  through  each  point  (xq,  Eq) , 
where  the  initial  wave  energy,  E , is  associated  with  x 
and  9q  on  the  initial  curve,  C.  Using  the  Runge-Kutta 
numerical  method,  Equation  (25)  is  integrated  with  the  same 
set  of  initial  values  of  n previously  defined  by  the  kinematic 
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consideration.  The  quanities  of  9,  99/3x^,  39/9x2/  |c|, 

3|c|/3x^  and  3|c|/3x2  are  determined  from  the  solutions  of 
Equation  (15)  as  functions  of  x in  the  integration  of 
Equation  (25). 

C.  Directional  Wave  Spectra 

The  waves  are  assumed  to  be  driven  by  wind  with  velocity, 
W,  which  is,  shown  in  Figure  1,  in  the  direction  of  the 
x^-axis.  The  most  convenient  way  of  describing  directional 
waves  is  the  directional  wave  spectra.  In  this  particular 
study  the  directional  wave  number  spectra,  ¥((£(,9),  is  used 
in  which  wave  energy  is  expressed  in  terms  of  a wave  whose 
local  wave  number  is  [ ic  | and  local  direction  of  propagation 
with  respect  to  wind  velocity  is  the  angle  of  9 . Up  until 
now  the  theoretical  derivation  for  the  directional  wave 
spectra  is  still  not  available  because  the  mechanism  involved 
in  the  generation  of  waves  by  wind  is  not  well  understood. 
Nevertheless,  for  practical  reasons,  there  are  some  semi- 
empirical  und  semi-theoretical  directional  spectral  forms 
which  have  been  proposed  and  studied.  The  majority  of  them 
are  the  mathematical  models  which  describe  waves  for  the 
region  from  which  data  are  collected.  It  is  not  the  purpose 
here  to  discuss  these  directional  wave  spectra.  For  more 
information  one  may  refer  to  Chen  (1972) . 

The  general  form  which  K(]kj,6)  is  assumed  is  based 
upon  the  assumption  of  separation  of  variables,  i.e, 

'M|k|,9)  = iM|)c|)  A ( | k | , 9 ) (32) 

where  ( | k | ) is  the  integrated  one-dimensional  wave  number 
spectra  and  A(|k|,9)  is  the  angular  spreading  function.  By 
satisfying  the  dimensional  analyses  of  Kitaigoroaskii  (1961) 
and  the  equilibrium  range  argument  of  Phillips  (1958),  t(|k  ) 
has  the  form  of 


where  the  value  of  q is  given  by  Equation  (36)  . By  the 
definition  of  wave  energy  density,  the  initial  wave  energy, 

— V t 

Eq,  for  the  wave  with  wave  number  vector,  kQ,  in  the  case 
with  no  current  can  be  obtained  by  using  Equations  (36)  and 
(37)  as 


VV  = V^ol 'V  - * <l*ol'0)  l*ol 

0.001  r(q  + 1) 

/?  r (q  + 1/2)  |ko|3 

• e-°-74  (9/lw|2|^0l)2  cos2’  ^ 


where 

= 1Q-0 . 3788  S W 1 (g/|kQ|)  1/2  + 1.059  (39) 

Equations  (38)  and  (39)  are  used  to  evaluate  Eo(kQ)  for  the 
numerical  calculations. 

IV.  NUMERICAL  RESULTS 

Contour  wave  energy  plots  of  the  directional  wave 

spectra  with  k,  or  |k|  and  9,  at  x = (50  km,  5 km)  are  made 

for  one  case  of  current  gradients.  Namely,  by  referring  to 

Equations  (22)  and  (23)  , current  gradients  are  chosen  to 
-5  -5 

have  a = 2 • 10  and  b = 2 • 10  . The  results  produced  by 

the  interactions  between  steady  non-uniform  currents  with 
these  current  gradients  and  wind-generated  gravity  waves  are 
given  in  Figure  2.  A reference  case  of  no  current  is  also 
given,  but  in  Figure  3.  The  active  wind  speed,  |w|,  in  all 
of  these  cases  is  chosen  to  be  10  meters  per  second.  In 
order  to  appreciate  the  results,  one  has  to  compare  the 
displacements  of  the  energy  contours  from  one  case  to 
another.  For  this  purpose  several  energy  contours  of  the 
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I 


same  set  of  values  are  plotted  on  each  figure.  The  results 
for  these  plots  are  obtained  by  utilizing  the  CDC  1604 
computer . 


J 


The  most  direct  computational  method  to  obtain  results 

for  the  contour  plots  is  the  use  of  root  subroutines,  as 

■+ 

shown  in  the  Appendix,  to  calculate  E = E(x,  k,  xQ,  kQ) 
where  | icQ  | and  n in  xQ  = (o,n)  are  the  independent  variables 
and  E is  the  solution  in  Equation  (25)  with  0,  39/3x^, 

30/3x2 , |c|,  3|c|/3x^,  and  3|c|/3x2  determined  by  the 
solution  of  Equation  (15),  8 = 9 (x,  xQ,  9 ) . For  each  value 
of  0Q,  the  value  of  xQ  is  varied  in  order  that  the  wave 
energy  characteristics  of  Equation  (25)  passes  through  the 
selected  space  point,  x,  and  the  value,  | icQ  | , is  also  varied 
so  that  wave  energy,  E,  at  this  space  point,  x,  lies  on  the 
wave  energy  contour;  for  each  space  point,  x,  along  the 
characteristic  curve  of  Equation  (25)  n is  varied  so  that 
the  characteristic  curve  of  Equation  (15)  passes  through  the 
space  point,  x,  of  the  integration  for  Equation  (25) . 

Furthermore,  to  decrease  the  time  required  on  the 
computer,  a generalized  Newton's  method  of  approximation  is 
used  in  addition  to  root  subroutines  to  obtain  solutions.  A 
listing  of  the  Fortran  program  is  given  in  the  Appendix. 

The  main  routine  PDEQ  integrates  Equation  (25)  for  initial 
values  of  n and  values  of  | ic  | ; the  subroutine  DISP  integrates 
Equation  (15),  calculating  9,  39/3x^,  30/3x2,  |c|,  3|c|/3x^ 
and  3|c|/3x2  at  each  space  point,  x,  in  the  integration  of 
Equation  (25)  by  the  main  routine. 

V.  DISCUSSION 

In  the  processes  of  solving  Equations  (15)  and  (25), 
which  are  the  first  order  partial  differential  equations,  by 
the  method  of  characteristics  numerically,  one  notices  that 
the  current,  U,  is  required  to  be  specified  initially,  on 
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the  initial  curve,  C,  as  well  as  along  each  characteristic 
curve.  Yet  this  specified  current  field  has  to  satisfy  the 
condition  of  continuity.  Equation  (9).  Because  of  these  two 
conditions  which  are  not  the  constraints  of  the  method,  the 
first  way  to  deal  with  the  problem  will  require  the  compu- 
tation of  wave  field  for  all  specifiable  current  fields  in 
the  ocean  which  satisfy  Equation  (9)  under  various  steady 
active  wind  conditions.  This  is  quite  a task;  even  the 
initial  reference  wave  field  is  specified.  The  second  way 
to  deal  with  the  problems  will  be  to  specify  the  wave  field 
and  to  solve  for  current  field  inversely  under  various 
steady  active  wind  conditions.  However,  the  processes  of 
solving  Equations  (15)  and  (25)  by  this  method  of  character- 
istics will  require  k(x)  and  wave  energy,  E(k),  to  be 
specified  initially  as  well  as  along  each  characteristic 
curve.  The  most  convenient  way  of  satisfying  these  require- 
ments is  by  using  the  remote  sensing  techniques  which  can 
provide  information  over  a large  area  synoptically . 

Nevertheless,  a few  cautions  are  warranted  here. 

First,  through  the  analysis  the  linear  dispersion  relation- 
ship as  shown  in  Equation  (6)  has  been  used.  This  implies 
that  the  weak  nonlinear  dynamic  interaction  is  neglected. 
However,  recent  theoretical  analyses  by  Huang  and  Tung 
(1976a  and  1976b)  show  that  the  reported  10%  increase  in 
wave  number,  | ic  | , as  compared  with  that  derived  from 
Equation  (6)  by  Longuet-Higgins , et  al . (1963)  can  be 
produced  by  the  directional  nonlinear  wave-wave  interactions. 
Improvements  can,  of  course,  be  made  to  include  such 
nonlinear  mechanisms  in  the  future. 

Second,  it  is  important  that  within  the  length  scale  of 
wave  decay  there  should  be  appreciable  changes  in  current 
velocity,  U.  Because  of  this,  capillary  waves  are  not 
considered  in  this  study.  Yet  if  we  concentrate  on  gravity 
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wave  range,  the  minimum  wave  length  of  0.3  meters  will  have 
a decay  time  of  2.486  • 10 J seconds  and  decay  distance  of 
1700  meters  for  kinematic  viscosity  of  9.17  • 10-7  seconds 
per  meter  square. 

Finally,  the  energy  balance  equation,  as  shown  by 
Equation  (11) , is  obtained  by  integrating  Navier-Stokes 
equation  over  the  depth  of  water  without  considering  the 
problem  of  wave  breaking.  This  wave  breaking  phenomenon 
will  cause  serious  mathematical  difficulties  in  the  case  of 
a saturated  ocean  encountering  opposing  currents  where  the 
vigorous  wave  breaking  will  drain  an  analytically  unaccount- 
able large  amount  of  energy.  Extra  caution  should  be  taken 
in  dealing  with  this  situation. 

VI.  CONCLUSION 

The  changes  in  wave  characteristics,  in  the  form  of 
directional  wave  spectra,  have  been  numerically  calculated 
as  the  directional  wind-generated  gravity  waves  encounter 
the  steady  non-uniform  two-dimensional  currents.  Using  the 
remote  sensing  techniques  which  are  capable  of  observing 
these  changes  in  waves  over  a large  area,  one  will  be  able 
to  infer  the  magnitude  and  the  direction  of  the  current 
encountered.  This  particular  method  of  inferring  current 
field  has  great  potential  benefits  in  global  applications. 
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